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Abstract 

We use molecular dynamics simulations to study the melting of gold icosahedral 
clusters of a few thousand atoms. We pay particular attention to the behavior of 
surface atoms, and to the equilibrium shape of the cluster. We find that the surface 
of the cluster does not pre-melt, but rather remains ordered up to the melting T^. 
However the increasing mobility of vertex and edge atoms significantly soften the 
surface structure, leading to inter- and intra-layer diffusion, and shrinking of the 
average facet size, so that the average shape of the cluster is nearly spherical at 
melting. 
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Gold particles consisting of tens to thousands of atoms have unique optical 
and mechanical properties and hold great promise as building blocks for nano- 
bioelectronic devices [1,2], catalysts [3], and sensors [4]. It is therefore natural 
that the physics and chemistry of these materials are a current research subject 
of great interest [5]. For future applications knowledge of the structure and 
stability of gold nanoparticles of different size and morphology is particularly 
important. 

While bulk gold has an fee crystal structure, the competition between bulk 
and surface energies in nanometer sized gold crystallites can result in several 
different competing structures [6,7]. Depending on cluster size and external 
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conditions transitions between these structures have been observed [8,7]. One 
such structure which has been observed both in simulations [9,10] and in ex- 
periments [11,12], is the "Mackay icosahcdron" [13,14], consisting of 20 shghtly 
distorted fee tetrahedra, with four {111} faces each, meeting at the center 
to form an icosahedral shaped cluster. The internal faces of the tetrahedra 
meet at strain inducing twin grain boundaries with hep structure, leaving the 
cluster with 20 external {111} facets. Theoretical models [15,16,17,18] have 
predicted different limits for the stability of such icosahedral clusters, and it is 
unclear whether their formation is an equilibrium or rather a kinetic process 
[11,17,18,19,20,21]. Nevertheless, it is natural to suppose that formation of 
this structure is related to the very high stability of the {111} external sur- 
faces. Simulations [22] and experiments [23] on bulk slab-like geometries with 
exposed {111} surfaces have shown that, unlike the {100} and {110} surfaces 
which melt below the bulk melting temperature T^, the {111} surface neither 
melts nor roughens but remains ordered up to and above Tm, and can in fact 
lead to superheating of the solid [25]. In light of this observation it is interest- 
ing to consider how the high stability of the {111} facets effects the melting 
and equilibrium shape of such icosahedral nanoclusters. 

In order to address this issue, we have performed detailed numerical simu- 
lations of icosahedral gold nanoclusters of a few thousand atoms, obtained 
by cooling from the melt. We pay particular attention to the behavior of 
the surface atoms and to the equilibrium shape. We find a sharp first order 
melting transition T^. Unlike earlier results on smaller cuboctahedral clusters 
[24], which include non {111} facets that pre-melt below Tm, we find no sur- 
face pre- melting of the {111} facets of our icosahedral cluster. Nevertheless, 
we find that there is a considerable softening of the cluster surface roughly 
~ 200 K below Ti„ due to the motion of atoms along the vertices and edges 
of the cluster. In this region we find both intra-layer and inter-layer diffusion 
of atoms, which increases considerably as is approached. The equilibrium 
shape progresses from fully faceted, to faceted with rounded edges, to nearly 
spherical just below Tm. Throughout this region, the interior atoms of the 
cluster remain essentially perfectly ordered, until Tm is reached. 

Using the many-body "glue" potential [26] to model interactions among gold 
atoms, we carry out molecular dynamics simulations, integrating the classical 
equations of motion with the velocity Vcrlct algorithm [27] with a time step of 
4.3 fs. The results presented here arc for a 2624 atom cluster, with a diameter 
of ~ 20 A, but we have also considered other sizes. We start our simulations at 
a high T — 1500 K > Tm, and cool using the Andersen thermostat method [27] 
to 1000 K. We then cool, in intervals of 100 K, down to 200 K, using 5 x 10^ 
steps (21 ns) at each temperature. Even though our N = 2624 atoms is not a 
"magic number" for a perfect icosahedral structure (the nearest such number 
being 2868), the cluster structure we find with this method is nevertheless 
clearly a Mackay icosahedron consisting of shghtly distorted tetrahedra with 
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Fig. 1. Potential energy vs. T for a 2624 atom gold cluster. The sharp jump at 
T = 1075 K indicates the first order melting transition. 

different number of atoms, but with a missing central atom. Similar icosahedral 
structures were obtained in runs with different particle numbers. Such a central 
vacancy, postulated for copper and aluminum but not for gold [28], has been 
considered previously [29] as a means of partially relieving the strain caused 
by the hep twin grain boundaries between the fee tetrahedra. 

To study melting and the equilibrium shape, we next heat the cluster up 
using constant temperature molecular dynamics[30]. To compute equilibrium 
properties, we take at each temperature 10^ steps (4.3 ns) for equihbration, 
followed by 10^ steps (43 ns) to compute averages. We take fine temperature 
increments in the vicinity of the cluster melting transition. In Fig. 1 we show 
our results for the average potential energy vs. temperature. We found the 
cluster to melt at = 1075 K, with a discontinuous jump in potential energy. 
Over the whole temperature range from 200 K to 1200 K gold atoms were never 
observed to evaporate from the cluster. 

To characterize the structure of the cluster, we measure the standard bond 
orientational order parameters Qq, Wq, Q4 and W4 [31] which are often used 
to distinguish between different phases of condensed materials [30,32]. These 
bond order parameters, designed to probe the degree and type of crystallinity, 
are sensitive to the orientational correlations of "bonds" , i.e. the vectors joning 
pairs of neighboring atoms. In the liquid phase, such correlations decay quickly 
with growing distance and the bond order parameters vanish. In crystalline 
solids, on the other hand, orientational bond correlations persist over large 
distances leading to order parameters with finite values. We refer the reader 
to the work of Ref. [31] for the definition of these quantities, and their values in 
common crystal structures. To distinguish between bulk and surface behavior, 
we first identify all atoms on the surface of the cluster, then the atoms in the 
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Fig. 2. Bond orientational order parameters, Qq,Wq,Q4,W4, vs. T, averaged over 
(a) bonds between interior atoms only, and (b) bonds between surface atoms only. 
For W4 only a few representative error bars are shown for the sake of clarity. 



first sub layer below the surface, and so on; the cluster has 9 such layers. We 
label as "interior" atoms those lying below the fourth sub layer (our results 
are essentially unchanged if we define the "interior" as all atoms below the 
first sub layer) . In Fig. 2a we plot the bond order parameters averaged over 
only bonds between the interior atoms. We see Q4, VK^ ~ at all T, while Qq 
and Wq take finite values for T < appropriate to the Mackay icosahedron. 
Qq and Wq remain essentially constant, decreasing only slightly just below Tm, 
indicating that the bulk ordering remain stable up until melting. In contrast. 
Fig. 2b shows the bond order parameters averaged over only bonds between 
the surface atoms. Again Q4, W4 ~ at all T, while Qq and Wq take finite 
values for T < Tm. However here we see a much more pronounced decrease, 
particularly in Qg starting well below Tm, until both vanish at melting (the 
finite value of Qq above melting is a finite size effect that decreases as the 
cluster size increases). Thus, while the surface remains ordered below T^, i.e. 
IQel; l^^el > 0, the surface order softens to a much greater extent than does 
the bulk as Tm is approached. 

Next we consider the diffusion of the surface atoms. In Fig. 3a we plot, for sev- 
eral different temperatures, the average mean square displacement (|Ar(t)p) = 
(|r(t) — r(O)p) vs. time t, where the average is over all the atoms which were 
initially on the surface of the cluster. At T = 600 K, diffusion is very slow, 
with displacements after 20 ns remaining less than one atomic separation. At 
T = 900 K, almost 200 K below Tm, diffusion is significant. At T = 1060 K, 
15 K below Tm, the mean square displacement saturates at large t, indicating 
that atoms now diffuse the entire length of the cluster. Fitting the early time 
linear part of these curves, (|Ar(t)p) ~ 6Dt, we plot the diffusion constant in 
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Fig. 3. (a) Mean square displacement for atoms initially on the surface, at various 
temperatures; (b) Diffusion constant D vs. T for atoms initially on the surface. 
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Fig. 4. For the indicated temperatures, top row: average cluster shape for a simu- 
lation time of 43 ns; middle row: histogram of the maximal local curvatures of the 
average shape; bottom row: ellipsoids indicate root mean square displacements of 
atoms on the cluster surface over a simulation time of 1.075 ns. 



Fig. 36. 
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The question thus arises how to reconcile this observed surface diffusion below 
Tin with the absence of surface melting that is indicated by the finite bond 
orientation order parameters of Fig. 2b. One possibility is that, as the temper- 
ature approaches T^, all surface atoms become more mobile but translational 
order is maintained due to the presence of a periodic substrate formed by the 
ordered sub-layers below the surface. Our simulations, however, point to a dif- 
ferent picture (see last row of figures shown in Fig. 4). For each atom initially 
on the surface of the cluster, we compute its average position f = (r), and 
its average displacement correlation matrix dap = {{r — f)a{r — f)^), where 
a, /3 — x,y, z and the (. . .) stand for averages over 25 configurations, sampled 
every 43 ps of simulation time. Taking the eigenvectors of dap and the square 
root of their corresponding eigenvalues to define the axes and principal radii 
of an ellipsoid, gives a convenient representation for the root mean square 
displacement of the atom. In the last row of Fig. 4 we plot these ellipsoids 
for each atom initially on the surface, centering the ellipsoid at the average 
position of the atom f . We show such plots for the same temperatures as in 
Fig. 3a. We clearly see that for T = 600 K and 900 K, the biggest ellipsoids 
are at the vertices and edges, while those for atoms in the facets are in general 
smaller. If we let the time that we average over increase, we find that the ellip- 
soids at the vertices and edges grow in size, corresponding to diffusion, while 
those in the middle of the facets stay approximately the same, corresponding 
to thermal vibrations without diffusion. A more detailed analysis shows that 
significant interlayer exchange takes place between the two topmost layers as 
much as 200 K below T^. In Fig. 4 one can see ellipsoids oriented normal to 
the cluster surface, corresponding to this interlayer diffusion. For the higher 
temperatures. T = 1060 K just below melting, and T = 1100 K just above 
melting, diffusion is pronounced throughout the entire surface. 

Finally, we consider the effects of the vertex and edge atom diffusion on the 
equilibrium shape of the cluster. Because of the small cluster size, the in- 
stantaneous shapes fluctuate significantly at high temperatures. But since our 
simulation algorithm conserves angular momentum, and the angular momen- 
tum is zero, our sample does not rotate as a whole. Hence we can compute a 
well defined average equilibrium shape at each temperature. We measure this 
equilibrium shape by averaging over the instantaneous shapes as follows. We 
divide space up into 842 approximately equal solid angles [33], corresponding 
roughly to the number of surface atoms. We then average the position of the 
surface atoms found in each solid angle over 1000 configurations sampled at 
equal times throughout the simulation of total time 43 ns. This defines the 
average radial position of the cluster within each solid angle, and hence the 
average cluster shape. In the top row of Fig. 4 we show the resulting equilib- 
rium shapes for several temperatures. We see that the shape is rounding out 
as the temperature increases, assuming a nearly perfect spherical shape above 
T 
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To quantify this, we compute the curvature distribution of the surface as 
follows. Using the average position of the surface in a given solid angle and its 
nearest neighbors, we fit to determine the best tangent plane to these points. 
Defining the normal to this plane as the z axis, we then find the best paraboloid 
that fits through the points. The principal curvatures of this paraboloid then 
give our approximation for the two principal curvatures of the surface at the 
given solid angle. We define k to be the maximum of these two principal 
curvatures. In the middle row of Fig. 4 we plot histograms of k as one varies 
over all the sohd angles defining the average surface. At T = 600 K and 900 
K we see a sharp peak at k = corresponding to the points on fiat facets, and 
a high K tail corresponding to higher curvatures on the edges and vertices. 
At T = 1060 K, just below Tm = 1075 K, the peak at k = has essentially 
vanished, and one has a broad distribution centered about ~ XjR with 
R — 21.5 A, the radius of the spherical liquid drop above Tm. At T = 1100 K, 
just above Tm, the distribution becomes very sharply peaked about k — 1/R. 

To conclude, we have found that gold nanoparticles of a few thousand atoms 
form a Mackay icosahedral structure, with missing central atom, when cooled 
from a liquid. Upon slow heating, we find that this bulk structure remains 
stable up to a sharp first order melting. The surface remains ordered with no 
pre-melting below Tm, however it softens considerably with increasing diffusion 
due to mobile vertex and edge atoms. As Tm is approached, this diffusion 
of edge atoms leads to significant shrinkage of the {111} facet sizes in the 
average cluster shape, leading to an almost spherical shape just below T^. 
In addition to the cluster of 2624 atoms reported upon here, we have also 
considered clusters of different sizes with 603 and 1409 particles. While the 
melting temperature was observed to increase with increasing cluster size, 
we continued to find the same general features, with the surface softening 
tracking the increase in T^. It would be interesting to know how this surface 
softening is related to morphological transitions observed in gold nanorods at 
temperatures below the melting temperature [30,34]. 

This work was funded in part by DOE grant DE-FG02-89ER14017. 
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